Cumulant ratios and their scaling functions 
for Ising systems in a strip geometry 
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We calculate the fourth-order cumulant ratio (proposed by Binder) for the two-dimensional Ising 
model in a strip geometry L x 00. The Density Matrix Renormalization Group method enables us 
to consider typical open boundary conditions up to L = 200. Universal scaling functions of the 
cumulant ratio are determined for strips with parallel as well as opposing surface fields. 
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Introduction. An universality principle is a corner- 
stone of contemporary theory of phase transitions. Ac- 
cording to this principle, the following sorts of quantities 
are universal: critical exponents, certain amplitude ratios 
and scaling functions ||I|. They differ each from other in 
their status. The (bulk) critical exponents are indepen- 
dent on boundary conditions, whereas two other groups 
are dependent. 

The critical exponents are known for many models 
(both exactly and approximately). The collection of re- 
sults available for amplitude ratios is also rich, but sig- 
nificantly smaller than for exponents; see in Ref. [Q for 
exhaustive information. 

Among amplitude ratios, so called cumulant ratios are 
of great importance. They supply some information on 
scaling functions (cumulants are proportional to deriva- 
tives of these functions at zero values of argument (s)); 
they measure deviation of magnetization fluctuations at 
criticality from gaussian distribution; moreover, they are 
closely related to some versions of renormalization group 
(it gains also reflection in terminology: cumulant ra- 
tios are customarily termed "renormalized coupling con- 
stant" in field theory). Cumulant ratios have also been 
used to locate the critical points and critical lines in many 
models §-||. 

Most results for cumulants have been obtained for 
(partially or completely) periodic boundary conditions. 
For most extensively studied the two-dimensional Ising 
model, numerous results are available ||2|,|5|-[l3|. However, 
there are also other, very natural "open" boundary con- 
ditions: "free" (no surface fields) , "wall-|— 1-" (infinite par- 
allel surface fields), "wallH — " (infinite opposing surface 
fields). For these open boundary conditions, the num- 
ber of results is very small. Only papers we know where 
such results are available are |g,Q|. Motivated by this 
situation, we state the aim of this paper: 

Calculation of universal cumulant ratios for the two- 
dimensional Ising model in a strip geometry under 
the following boundary conditions: "free", "wall-\ — \-", 



"wall+- ". 

We have calculated cumulant ratios using method 
called the Density Matrix Renormalization Group 
(DMRG). Since the DMRG is most powerful for open 
boundary conditions, it is particularly suited for our 
goals. 

Definition of cumulants. We consider the two- 
dimensional Ising system on a square lattice in a strip 
geometry {L is width of the strip and N is its length) 
with the Hamiltonian 
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where the first sum runs over all nearest-neighbour pairs 
of sites while the last two sums run over the first and the 
i-th column, respectively. H is the bulk magnetic field, 
whereas Hi and H]^ are the surface fields. H , Hi and 
-ffi are dimensionless quantities (all of them are mea- 
sured in units of J) . In the course of calculations of the 
cumulant in the termodynamic limit, two limiting pro- 
cesses are taken: T ^ Tc and L — > 00. In general, a 
value of the cumulant do depend on ordering of these 
limits 111] . In our paper we analyze so-called " massless" 
case (analogously as in |l|,|l|l): T = T^ ^ 2/ln(l + ^2) 
~ 2.269185 followed by L — > 00. Therefore, we do not no- 
tice the temperature dependence below. We also drop (as 
unnecessary) explicit dependence on surface fields until 
discussion of scaling functions. 

We consider the ratio of moments of magnetization 
proposed by Binder J2|. Definitions of cumulant ratios 
for a system in a strip geometry have been widely pre- 
sented in the literature P,p|,p|,pd]| . Let us first define 
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where M = ^^ Si is the total (extensive) magnetization. 
Then, the cumulant ratio Ajj in question is 
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lim L-'^Ut 
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An equivalent (but more convenient for us) formula 
for the above cumulant is as follows: Let A(L; H) be the 
largest eigenvalue of transfer matrix for the strip of width 
L {—T log X{L; H) is the free energy for one column of 
spins). Define 
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Then our cumulant is equal to pi 

Au = lim r{L) = lim mi{L)/2,Lml{L). (6) 

L — >oo L — >oo 

Our method of calculation is based directly on the def- 
inition (|^). We first find logarithms of the largest eigen- 
value \{L\H) for some values of H (at fixed L). Next, 
we calculate numerically derivatives (Q) and (||) , then the 
ratio r{L), and finally perform the extrapolation L — > oo. 

Some technical details of calculations. We use 
the DMRG method for calculations of log A(L; H). Orig- 
inally, this method has been proposed by White ||l^ for 
finding accurate approximations to the ground state and 
the low-lying excited states of quantum chains. Its heart 
is recursive construction of the effective Hamiltonian of 
a very large system using a truncated basis set, start- 
ing from an exact solution for small systems. Later, the 
DMRG was adapted by Nishino ||l6| for two-dimensional 
classical systems. The DMRG has been applied suc- 
cessfully to many different problems and now it can be 
treated as a standard method, which is very flexible, rela- 
tively easy to implement and very precise. For a compre- 
hensive review of a background, achievements and limi- 
tations of DMRG, see in Refs. ^^. 

A factor crucial for precision of DMRG is so-called 
number of states kept m, describing the dimensionality 
of effective transfer matrix |^,0 ; the larger number of 
states kept, the more accurate the value of the free en- 
ergy. Using m = 50 we can calculate the free energy 
with accuracy of the order 10~^^ for strips of width of 
the order L = 200. This is an one order more than size 
of systems available by exact diagonalization of transfer 
matrix. This fact is crucial for us, because of using the ex- 
trapolation procedures. In our calculations we apply the 
finite system, algorithm, developed by White for studying 
finite systems |l3] . Additional factor determining the ac- 
curacy of the method is number of sweeps, i.e. number 
of iterations made in order to obtain self-consistency of 
results. Our numerical experience shows that in most 
cases, it is sufficient to apply only one sweep (although 
in the "wallH — " case two sweeps are necessary - see be- 
low). 



In our calculations of cumulant ratios, we have also 
a factor limiting accuracy that is independent on the 
DMRG method: accuracy of numerical differentiation. 
In the procedure of numerical differentiation, a suitable 
choice of increment AH of an argument is of crucial im- 
portance. It is clear that AH should be taken as small 
as possible; on the other hand, due to finite accuracy 
in calculation of A, an error of difference quotient in- 
creases with decreasing AH. The increments used in 
our calculations have been determined as a compromise 
between above two tendencies. Additional factor deter- 
mining the accuracy of numerical differentiation, is a 
number of points used to calculate the derivative. We 
use formulas where a derivative is determined from the 
second-order Taylor expansion (i.e. we need n + 3 values 
of function for n-th derivative; this way, an accuracy is 
of the order 0{{AH)'^).) Therefore, the m2 was deter- 
mined from 5 points (3 points in symmetrical case, i.e. 
f{H) = f{—H)) and from 7 points for 7714 (4 points when 
the symmetry was present). 

We have tested correctness of our calculations in sev- 
eral ways. One of them was L-dependence on derivatives 
TO2 and rrn. Finite-Size Scahng (FSS) theory |l|,|l9[ pre- 
dicts the following dependence on the n-th derivative of 
the free energy as a function of the system size L: 
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where we have A = 15/8 and v = \ for the two- 
dimensional Ising model, d is a dimension of system in 
"finite-size direction" , i.e. it is a number of linearly in- 
dependent directions along which a size of the system is 
finite. For finite systems (for instance a torus), d is equal 
to space dimension of the system. In our case, the system 
is infinite in one direction (along the strip) and finite in 
the second direction (across the strip), so we have to take 
d = \. This assumption gives the following predictions 
for derivatives: 
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where p2 = 11/4 and p^ — 13/2. 

An extrapolation procedure has been performed with 
use of the powerful BST method |2^. 

Results: the "free" case. It corresponds to zero 
surface fields Hi = Hl = in the formula (|]). We have 
performed calculations for L in the range 160 < L < 200 
with step 10; these values of L were taken in all situ- 
ations. We took an increment of "bulk" magnetic field 
AH = 5 X 10~^, m — 50 and one sweep. The results are 
listed in the Table. 

As a byproduct, we have tested the FSS predictions for 
L-dependence of derivatives m2 and TO4. Values of cor- 
responding exponents (see Eq. (||) ) are: p2 = 2.7495(3) 
and p4 = 6.50(3), so predictions of FSS are confirmed in 
excellent manner. The same conclusion is true in next 
two situations. 



As another test of correctness (and quality) of DMRG 
results, we have calculated ratios by immediate numer- 
ical diagonalization of transfer matrix for 10 < L < 18 
(L even; these values of L have also been used in the 
next cases). We proceeded as above, i.e. by calcula- 
tion of logarithm of the largest eigenvalue for some val- 
ues of bulk field ff, followed by numerical differentia- 
tion of f{H) and computation of ratio and extrapola- 
tion, without any "renormalization" . We took an incre- 
ment AiJ = 10"^. We have obtained Au = -1.094(1); 
P2 = 2.746(1); p4 = 6.46(1). It is seen that the results 
are fully consistent with the DMRG calculations but less 
precise; we have the same situation for two other bound- 
ary conditions. 

Results: the "wall++" case. The "wall++" 

boundary condition corresponds to the assumption that 
all boundary spins have the same value and sign. It is 
equivalent to putting Hi = H^ = oo in (|l|). Numerical 
experience suggests that it is sufficient to take Hi — 10 
- for larger values of Hi the changes of the free energy 
are negligible [£1[ . 

The "wall-|--|-" configuration is more intricate, from nu- 
merical point of view, than "free" system. The compli- 
cation is due to the fact that, for parallel surface fields 
of the same sign, the maximum of the free energy f{H) 
does not appear for H = but it is shifted to a cer- 
tain non-zero value Hq(L). This phenomenon is called 
the capillary condensation pl| , |2^ . In order to calculate 
derivatives and ratios at zero magnetization (i.e. at the 
maximum of the free energy), we first have to find its 
position Ho{L). FSS predicts |22| the following depen- 
dence: Hq{L) ^ L^^l^ . From our DMRG calculations, 
we have obtained the value ^jv = 1.8749(2). 

For the "wall-|--|-" configuration the free energy is 
not longer a symmetric function of the bulk field _ff, 
so we have been forced to calculate m2 from 5 points 
and ?7i4 from 7 points. We have taken an increment 
AiJ = 5 X 10^®, m = 50 and one sweep. The results 
are presented in the Table. For exponents of m2 and 7714, 
we have obtained: pa = 2.7504(3) and pi = 6.5024(3). 
The precision of these results is a little bit less than in 
the "free" case, although still very satisfactory (three sig- 
nificant digits). However, it should be stressed that here 
we must do much more numerical computations than in 
the "free" case, so some lack of precision is inevitable. 

Exact diagonalization of transfer matrix gave the fol- 
lowing values: Au = 0.45(4); p2 = 2.75(1); pi = 6.5(2). 

Results: the "wallH — " case. One of important 
physical implications of "H — " boundary condition is the 
presence of an interface between "-I-" and "— " phases 
in the system. It causes large fluctuations, which have 
an implication in numerical practice, namely, two sweeps 
are necessary to ensure self-consistency of results. In our 
calculations, the value of surface field Hi ~ 100 an in- 
crement AiJ = 2 • 10~^, and m = 40 were taken. The 
results are listed in the Table. Values of exponents are: 



P2 = 2.7502(2), /94 = 6.502(2). 

The procedure of exact diagonalization of transfer ma- 
trix gave the following values: Au = —0.305(2); p2 — 
2.755(1); /94 = 6.50(2). 

As a matter of some interest, let us remark that for the 
"wallH — " boundary condition the L dependence is much 
weaker than for the "free" and "wall-|-+" situations. 

Table. Values of cumulant ratios for some values of L 



L 


r{L), free 


r{L), wall+-K 


r(i), wallH — 


160 


-1.098525 


0.462556 


-0.304831 


170 


-1.098234 


0.462225 


-0.304859 


180 


-1.097964 


0.461859 


-0.304883 


190 


-1.097723 


0.461525 


-0.304902 


200 


-1.097481 


0.461133 


-0.304915 


00 


-1.0932(3) 


0.455(2) 


-0.3050(1) 



Scaling functions for ratios. The "wall -I- -f / H — "- 
type conditions can be treated as limiting case of the 
system with equal finite parallel/antiparallel (+ + / H — ) 
surface fields. Another limiting case is the "free" bound- 
ary condition, where the values of surface fields are set 
to zero. One can expect that for intermediate situa- 
tions, i.e. finite values of boundary field iJi, cumulants 
would be smooth functions of Hi . Particularly interest- 
ing are scaling properties of these functions. The scaling 
theory predicts that at criticality, the system depends 
only on one variable, namely dimensionless combination 
C = LHf |23]. In the other words, we expect that the 
cumulant r should depend on the surface field Hi and 
strip width L only by combination C. We have calculated 
r{L,Hi) for both "++" and "H — " boundary conditions 
for L — 40, 80, 120 using m = 40, and at full range of 
scaling variable. The results are presented on Figs, la 
and lb. It is seen that scaling properties are confirmed 
in excellent manner. 

Limiting values of these functions (i.e. for C = and 
^ ^ cxd) are fully consistent with our more precise calcu- 
lations, although convergence of ratio to its limit value 
is much faster for "wallH — " than for "wall-!--)-". As far 
we know, scaling functions for cumulants have been al- 
most not studied so far. Only exception is paper [E4|; 
however, the authors consider scaling functions different 
from ours. 

Summary. We have calculated cumulant ratios for 
Ising strips with three natural boundary conditions, al- 
most not studied so far: "free", "wall-|--|-" and "wallH — " 
situations. We have applied the Density Matrix Renor- 
malization Group method followed by numerical differ- 
entiation and extrapolation L —* 00. We claim that our 
results are very precise (three or four significant digits). 
The precision is comparable with three other "top qual- 
ity" methods used in similar calculations: Monte Carlo 
p2[ , some versions of Renormalization Group |l3[ and 
analysis of high-temperature series pOt] . 
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FIG. 1. The cumulant ratio as a function of scaling vari- 
able (^ = LH^: a) the "++" boundary condition; h) the "H — " 
boundary condition. Notice much weaker L dependence for 
"H — " than for "++" situation. Moreover, the convergence of 
ratio to its limit value, when (" ^ oo is much faster for "-|-— " 
than for " ++" . For the latter case, the saturation is achieved 
for C, ~ 500 (outside the range of plot). 

We have also calculated the quantity which apparently 
escaped the attention so far (at least for boundary con- 
ditions considered by us), namely, the scaling functions 
for cumulants. Such functions give information of how 
finite surface fields influence values of cumulants. This 
influence is significant - in one case ("-I--I-") even the sign 
of cumulant changes upon growth of surface field. We do 
not know how fundamental this phenomenon is. At first 
glance, it seems to be related to the lack of symmetry 
(i.e. f{Ho + H) ^ f{Ho — H)) and to a non-zero value 
of the third derivative of the free energy / with respect 
to H at Hq. However, the other explanations are not 
excluded and we will discuss it elsewhere. 

Natural lines of continuation of our investigations are: 
testing of universality of cumulants and scaling functions 
(for other models in the two-dimensional Ising univer- 
sality class, for example the hard squares model) and 
calculation of higher cumulants. This work is currently 
in progress. 

Most of our numerical calculations were performed on 
Pentium II machines under Linux. We used ARPACK 
package to calculate eigenvalues. 
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